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ABSTRACT 


A  method  is  presented  for  solving  hydrodynamic  prob¬ 
lems  involving  large  distortions  and  compressions  of  the 
fluid  in  several  space  dimensions.  The  calculation  procedure 
introduces  finite  difference  approximations  to  the  differ¬ 
ential  equations;  the  solution  in  practice  is  carried  out  by 
means  of  high-speed  electronic  computers.  The  paper  dis¬ 
cusses  a  number  of  characteristics  of  the  method  and  il¬ 
lustrates  these  by  presenting  results  of  representative  calcu¬ 
lations  . 


-  3  - 


CONTENTS 


Abstract 

Introduction 

I.  The  Method 

A.  General  Description 

B.  PIC  in  One  Space  Dimension 

C.  Discussion  of  1-D  PIC  Method 

II.  Extension  of  Method  —  Two  Space  Dimensions 

A.  The  Plane,  Cartesian,  Two-Dimensional  Box 

B.  Problems  in  Other  Coordinate  Systems 

in.  One-Dimensional  Tests 

A.  Simple  Steady  State  Shock 

B.  The  Closed  One-Dimensional  Box 

C.  Expanding  Spheres 

IV.  Two-Dimensional  Calculations 

A.  The  SUNBEAM  Code 

B.  The  KAREN  Code 

C .  Cylindrical  Coordinates 

Appendix  I.  Other  Methods  for  Solving  2-D  Problems 

Appendix  n.  An  Alternative  Derivation  of  the  PIC 
Difference  Equations 


Page 

3 

7 

8 
8 
8 

12 

23 

23 

27 

29 

29 

30 
32 

34 

34 

35 

35 

36 


39 


-  5  - 


INTRODUCTION 


Realistic  studies  of  the  complicated  dynamics  of  compressible  fluids 
are  being  made  possible  by  recent  development  of  large  high-speed  com¬ 
puting  machines.  Methods  are  well  known  for  treating  problems  dependent 
upon  one  space  coordinate  only  ("1-D"  problems).  For  those  dependent  upon 
two  or  more  space  coordinates  ("2-D"  problems,  etc.),  several  methods  of 
treatment  have  been  used  with  success  in  restricted  classes  of  problems. 
Features  of  these,  pertinent  to  this  report,  are  discussed  in  Appendix  I. 

It  is  our  purpose  here  to  present  and  discuss  a  method*  for  solving 
2-  and  3-D  problems  which,  in  a  number  of  preliminary  studies,  has  ap¬ 
peared  successful  for  solving  problems  involving  large  distortions  in  sys¬ 
tems  with  several  fluids. 

We  acknowledge  with  pleasure  the  stimulation,  encouragement,  and 
valuable  criticisms  that  have  come  from  numerous  discussions  with  Garrett 
Birkhoff,  Eleazer  Bromberg,  Rolf  Landshoff,  Robert  D.  Richtmyer,  A.  H. 
Taub,  and  Stanislaw  Ulam. 


*F.  H.  Harlow,  J.  Assoc.  Comp.  Mach.,  4,  137  (April,  1957). 
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Chapter  I 


THE  METHOD 


A.  General  Description 

The  calculation  procedure  which  is  here  described  combines  features 
of  several  of  those  discussed  in  Appendix  I.  The  space  occupied  by  the 
fluid  is  divided  into  a  network  of  fixed  cells  through  which  the  fluid  moves. 
The  fluid  within  these  cells  is  represented  by  particles,  each  of  which  carries 
a  fixed  mass  of  fluid.  Their  coordinates  vary  with  time  as  they  move 
through  the  mesh  of  cells.  Thus  we  have  a  Lagrangian  coordinate  system 
superimposed  upon  an  Eulerian  one. 

For  convenience,  we  refer  to  the  procedure  to  be  described,  as  the 
"Particle— in -Cell  Method,  "  abbreviated  "PIC . " 

B.  PIC  in  One  Space  Dimension 

Characteristics  of  method  are  most  easily  demonstrated  in  application 
to  1-D  problems,  but  PIC  attempts  competition  with  other  methods  only  for 
2-  and  3-D  calculations. 

Assume  only  one  material  to  be  present  in  a  closed  one -dimensional 
box.  The  system  of  equations  which  we  wish  to  solve,  subject  to  initial  and 
boundary  conditions,  is: 
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(3) 
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p  =  Up,  D  (4) 

where 

p  =  density 
u  =  velocity 
p  =  pressure 

I  =  specific  internal  energy 

1  2 

E  =  specific  energy  =  I  +  — u 

Relation  (4)  represents  the  equation  of  state  of  the  fluid.  The  absence  of 
heat  conduction,  viscosity,  or  other  dissipative  mechanism  means  that  the 
true  solution  may  possess  or  develop  mathematical  singularities.  This 
possibility  is  ignored  at  first. 

The  box  is  divided  into  cells  of  equal  length,  s,  labeled  with  index 
c  (c  =  1,  2,  .  .  .,  L).  Into  this,  particles  of  equal  mass,  m,  are  arranged 
with  number  Nc  in  each  cell  in  such  a  manner  that  the  initial  density  profile 
is  represented. 

The  calculation  to  advance  the  configuration  in  time  proceeds  in  steps 
or  cycles,  each  of  which  calculates  the  desired  quantities  for  time  t  +  tit  In 
terms  of  those  at  time  t  (an  ’’explicit"  time  advancement  procedure).  This 
introduces  nothing  novel,  however,  as  long  as  tit  is  small  enough  to  satisfy 
criteria  discussed  later.  Our  principal  attention  is  devoted  to  the  spatial 
features  of  the  calculation. 

A  division  of  each  cycle  of  calculation  into  three  phases  is  convenient. 

In  Phase  I  the  Eulerian  field  functions  are  changed,  neglecting  transport  due 
to  fluid  motion.  Phase  n  accomplishes  the  motion  of  the  particles  to  their 
new  positions.  The  transport  corrections  are  accomplished  in  Phase  in. 
Generally,  a  Phase  IV  is  also  added  in  which  various  diagnostic  functionals 
of  motion  are  calculated  for  that  cycle. 

Phase  I 

-  N  m 

Q 

Each  cell  has  density,  pc  =  — - — ,  specific  internal  energy,  I,,,  velocity, 
u^,  and  pressure,  pc  =  f(pc,  l£).  All  of  these  are  the  values  at  the  cell 
centers,  but  are  assumed  to  hold  for  a  particle  anywhere  in  a  cell.  Dropping 
the  transport  term  in  Equation  (2),  and  transforming  it  to  our  finite  difference 
form, 
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N  m  8u 
c  _ c 

s  8t 


2s  (Pc+1  Pc-l) 


A  similar  treatment  of  Equation  (3)  leads  to 


N  m  /  81  8u 

C  I  _ c  _ c 

s  \  8t +  Uc  8t 


-  ~ti  [  wM1  -  <p“)c.i] 


These  equations  can  be  combined  and  rearranged: 


8t  2N 


m  (Pc-1  Pc+l) 


=  ITS  [Pc-l(Vl  -  uc)  -  pc+l(uc+l  -  uo) 


When  the  intitial  and  boundary  conditions  have  been  specified,  then  the 
quantities  on  the  right  sides  are  all  known  at  time  t  and  the  time  deriva¬ 
tives  can  be  calculated. 

Let  t  =  ndt.  Then 


u  =  u  +  6t 
c  c 


i  -f  +  ot  hr 

c  c  \  8t  / 

where  the  superscript  represents  the  time  cycle  number  and  the  tilde  sym¬ 
bolizes  the  fa ct  that  these  advanced  time  quantities  are  yet  to  be  corrected 
by  the  transport  terms. 

An  alternate  form  for  the  energy  equation  can  be  obtained.  If  in 
Equation  (6)  we  make  the  change, 


9uo  1 8uc  1  r~  2  ,v 
'c  “it  <*0>  -<V 


and  shift  the  result  to  the  right  side  of  the  equation,  we  would  obtain,  in 
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place  of  Equation  (8), 


+  2N  m  ( Vl  Uc-1  Pc+1  Vl) 
c 


(8T) 


All  the  calculations  reported  in  this  paper  used  the  energy  equation  in  the 
form  of  Equation  (8).  Boundary  conditions,  the  same  for  either  form,  are 
discussed  later  in  connection  with  questions  of  conservation. 

Phase  n 

The  particles  are  moved.  Although  each  particle  in  cell  c  now  carries 
momentum  mUg,  it  is  moved  according  to  a  "velocity-weighting*'  procedure. 
The  effective  moving  velocity  is  obtained  by  a  linear  interpolation  according 
to  the  position  of  the  particle  between  the  centers  of  two  cells,  using  the 
velocity  at  each  center.  Labeling  the  particles  with  j,  we  write 


n+1  n  nm 

xj  “  x)  +  ueff  st  <10) 

The  value  of  u  „  is  computed  from  values  of  u11. 

eft  c 

Phase  HI 

The  transport  corrections  are  accomplished  in  this  phase.  For  those 
cells  which  change  particle  number  in  Phase  n,  additional  calculations  are 
required.  The  procedure  is  called  "repartitioning.  *'  Otherwise,  the  tilde 
values  of  Phase  I  become  final,  and  the  cycle  is  complete. 

Suppose  that  the  only  boundary  crossing  involves  a  particle  going  from 
cell  c  -  1  to  cell  c.  For  cell  c  -  1  the  only  change  recorded  is 

N11*!'  =  N*1  „  —  1.  Correspondingly  N„  is  increased  by  1.  Nc  ■  ^ c.  ^ 

c  —1  c  —1  '■* 

Furthermore,  the  velocity  of  cell  c  becomes 


n+1 


u 


N  u  +  u 
c  c  c— 1 

n 

N+1 

c 


(11) 
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-j  i  k  0  I 


(nof 


so  that  the  momentum  of  tiie  ^system  is  unchanged./  The  kinetic  energy  of 
the  system  drops,  however,  the  change  being 
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%hif 


The  specific  internal  energy  is  modified  by  the  transport  effect  and,  to  con¬ 
serve  total  energy,  is  increased  by  the  amount  of  loss  in  specific  kinetic 
energy. 


XI  rsw/ 

,»+l  .  V°  °-l,  l«  ke|  (is, 

C  -ft  *  !) 

When  a  cell  has  both  outgoing  and  incoming  particles,  the  latter  repartition 
only  with  those  which  remain  after  the  outgoing  ones  have  left.  Otherwise, 
an  unreal  transport  of  energy  is  added. 


C.  Discussion  of  1-D  PIC  Method 


1.  Energy  and  Momentum.  The  difference  equations  of  Phase  I  can  be 
obtained  in  several  ways  and  written  in  several  forms.  A  discussion  is 
given  in  Appendix  n.  In  particular,  Equation  (8)  is  more  complicated  'than 
the  analogous  equation  one  would  obtain  by  starting  with  the  differential 
equation 


I 

at 


(14) 


instead  of  Equation  (3).  The  reasons  for  our  choice  are  brought  out  by 
looking  at  the  problem  of  over-all  conservation  of  energy  in  Phase  I. 
The  total  energy  of  the  system  at  time  t  is 


(15) 


Likewise,  by  the  end  of  Phase  I  it  has  changed  to 


N  m 
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(16) 


-  12  - 


the  total  change  being 


f-£n  =r  6E 
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Define: 


Ncm(“c  -  Uo) 
c 

These,  with  Equations  (6)  and  (9),  become 

L 

6E  =  6D  +^|  2  [(pu)o_1  -  (PU)0+1J 


c 

>  "l 


+  V 


(18) 


/,K'fU3 

/pi >  o'  " e<) 


(19) 


2  /  8u  N 

c  v  ‘ 


(20) 


In  terms  of  two  fictitious  cells  beyond  the  system  (  c  =  0  and  c  -  L  +  1) 
Equation  (19)  can  be  written 


6E  =  6D  +  St 


(PU)Q  +  (pu)1  (pu)L  +  (pu)L+1 


(21) 


The  cancellation  of  terms  in  pairs,  exploited  here,  would  not  have  occurred 
had  we  started  from  Equation  (14)  in  the  difference  equation  derivations. 

The  important  quantity,  6D,  represents  the  only  discrepancy  in  energy 
of  the  entire  cycle  calculation,  since  Phases  n  and  in  conserve  energy 
rigorously.  This  nonconservation  would  not  have  appeared  had  Equation  (8*) 
been  used  fo ^Internal  energy  calculations.  Since  SD  0,  the  total  observed 
energy  rises  monotonicaliy  above  that  which  theoretically  is  present.  This 
fact  can,  indeed,  be  generalized  to  include  multidimensional  applications  of 
PIC  with  various  boundary  conditions,  conservative  or  not.  SD  is  an  im¬ 
portant  Phase  IV  functional  to  be  calculated  since  its  disagreement,  beyond 
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r£X'  7/  ' 
~  Sb  u 


x"  -J  //c  yn  £_  Ut  J 

machine  round-off,  with  the  observed  discrepancy  is  a  powerful  indicator  of 
almost  any  error,  mechanical  or  human. 

Knowledge  of  6D  is  of  value  also  for  the  purpose  of  limiting  the  size  of 
<5t.  Since  6D  is  proportional  to  (fit)  2,  its  cumulative  effect  over  a  fixed  time 
interval  is  proportional  to  fit.  The  time  interval  should  be  small  enough  so 
that  the  over-all  discrepancy  through  a  problem  is  insignificant.  Effectively, 
this  means  ufit/u  «  1  is  required,  since  the  ratio  of  the  specific  discrepancy 
to  the  specific  kinetic  energy  is  just  (ufit/u)2 .  Notice  that  when  the  velocity 
starts  from  zero,  then  in  the  first  cycle,  all  the  kinetic  energy  that  is  pro¬ 
duced  goes  into  discrepancy! 

Tests  have  shown  that  when  fit  is  small  enough  to  produce  negligible 
values  of  5D,  then  all  inaccuracies  due  to  time  differencing  are  negligible. 

In  this  case,  calculations  using  either  Equation  (8)  or  (8*)  are  equivalent,  and 
the  more  convenient  form  may  be  employed. 

Equation  (21)  shows  how  the  boundary  conditions  are  to  be  chosen  for 
conservation  of  energy:  the  average  value  of  pu  (the  energy  flux)  must  vanish 
at  each  boundary.  With  Pq  =  p^  and  uq  =  — u^  the  boundary  is  reflective. 

No  other  arrangement  is  reasonable  in  this  case. 

,  -  Since  Phases  n  and  HI  rigorously  conserve  momentum,  we  may  derive 
from  Equations  (5)  and  (9)  that  the  change  in  momentum,  6M,  in  a  cycle  is 


given  exactly  by 


6M  =  fit 


M  "A!1 


P0+Pl 


PL  +  *L+ 1 


(-%,»)( 

(22) 


Comparison  of  the  observed  change  with  this  theoretical  change  in  each  cycle 
makes  another  good  calculation  check,  since  the  agreement  should  be  to  within 
machine  round-off  errors. 

A  complication  arises  in  the  Phase  I  calculations  wherever  a  cell  be¬ 
comes  empty.  Equations  (7)  and  (8)  cannot  be  used  for  such  cells;  indeed, 
one  need  not  bother  to  calculate  in  them  at  all  in  Phase  I.  But  unless  special 
care  is  devoted  to  the  adjacent  cells,  the  system  will  suffer  significant 
energy  discrepancies.  A  reasonable  procedure  which  is  consistent  with  the 
y  conservation  laws  is  as  follows.  Let  cell  c  be  empty.  When  calculating  in 
/  cell  c  -  1,  consider  that  cell  c  has  pc  =  “Pc_j»  uc  =  +uc-l*  w^ien  ca^cu~ 

C  /  f  lating  in  cell  c  +  1,  consider  that  cell  c  has  pc  =  —  Pc+1»  uc  =  +uc+l*  Sim- 
>X>’ilar  procedures  hold  in  the  applications  of  PIC  to  2-D  problems. 

*>  -/Phi if  »><•<••>'  t>*  5/' •V'/V'!  {.■]■  V-'-'.-M/ 

2.  Diffusion  and  Accuracy.  Certain  other  important  characteristics  of 
PIC  are  shown  by  the  following  analysis. 

*  We  define  pn  as  the  average  mass  per  unit  distance  at  time  nfit  and 

c+a 
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position  (c  +  a)s.  Then  to  lowest  order  in  the  cell  size  and  zero  order  In 
<5t,  the  result  of  all  three  phases  of  the  calculation  can  be  summarized  by 
the  equations 


31 


at  2p  s 
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+  2p  8  PC~i  uc-|  (Vl  uc) 
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+  p  1  u  ,  , 

c~h  \  Pc S 


3  P, 


s  at  ^c-g  uc-g  Pc+J  uc+2 


(25) 


where  it  is  assumed  that  all  the  velocities  are  positive  in  the  direction  of 
increasing  cell  index.  Also, 


u 

c 


+  u 


c±l 


Equation  (25)  is  clearly  in  conservation  form, 
put  into  the  completely  conservative  forms 


(26) 

The  other  two  can  also  be 
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by  some  simple  algebraic  manipulations. 

Now  expand,  by  Taylor  series,  the  right  side  of  Equation  (23)  about  the 
center  of  cell  c.  Dropping  the  index  we  get 


(27) 
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Similarly  from  Equation  (24), 
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and  from  Equation  (25) 


BP  = 

at 


(28) 


(29) 


where 

X '  =  ~pus  (30) 

u 

From  these  results,  many  of  the  most  striking  characteristics  of  PIC 
can  be  learned.  Of  particular  interest  is  the  fact  that  repartitioning  does 
incorporate,  to  lowest  order,  the  correct  transport  terms  into  the  equations. 
The  next  order  term  In  the  momentum  equation  has  exactly  the  form  of  a 
"true"  viscosity  term  with  coefficient  of  viscosity  X1  =  £  Pus  •  The  corre¬ 
sponding  term,  X’  (0u/ax)2,  in  the  energy  equation  also  has  the  proper  form 
of  a  '♦true"  viscosity  effect. 

The  form  of  the  one-dimensional  effective  viscosity  X'  bears  a  striking 
resemblance  to  the  "linear  viscosity"  introduced  by  Landshoff*: 


where  c  is  the  local  sound  speed.  The  resemblance  to  the  classical  arti¬ 
ficial  viscosity  of  von  Neumann  and  Richtmyer**  is  less  close,  as  theirs 
possesses  a  velocity  gradient  factor. 

The  X'  terms  have  several  effects  on  the  calculation.  Most  important, 
they  symbolize  the  dissipation  which  is  present  in  regions  of  rapid  change, 
such  as  regions  which  would  contain  shocks  in  the  true  solution.  They  allow 
the  calculation  to  reproduce  correctly  the  macroscopic  features  of  a  shock 
without  requiring  special  care  concerning  Its  microscopic  features. 

The  X'  terms  also  symbolize  inaccuracy  in  the  calculations.  But  at  the 


*  R.  Landshoff,  Los  Alamos  Scientific  Laboratory  Report  LA-1930,  "A  Nu¬ 
merical  Method  for  Treating  Fluid  Flow  in  the  Presence  of  Shocks,  "  De¬ 
partment  of  Commerce,  Office  of  Technical  Services,  Washington  25,  D.C.  (1955). 

**j.  von  Neumann  and  R.  Richtmyer,  J.  Appl.  Phys.,  21,  232  (1950). 


same  time  they  suggest  a  cure,  since  the  effect  of  the  terms  can  be  made 
cumulatively  as  small  as  desired  by  decreasing  the  cell  size.  Various 
specific  criteria  can  be  established  of  the  following  sort.  Suppose  that  p  and 
u  are  everywhere  constant  in  space.  Then  the  condition  that  the  diffusion 
term  be  negligible  compared  with  the  transport  term  in  the  energy  equation 
is  that 


(31) 


The  right  side  is  a  measure  of  the  distance  over  which  the  slope  changes 
appreciably.  By  similar  conditions  applied  to  other  field  functions  one  may 
arrive  at  the  general  statement: 

PIC  will  produce  relative  errors  in  the  equations  of  motion  which  are 
roughly  equal  to  the  ratio  of  the  cell  size  to  the  smallest  physically  signif¬ 
icant  dimension  of  the  system. 

There  is  some~hope  that  this  statement  is  pessimistic.  Certainly  some 
functionals  of  the  motion  will  be  more  truly  represented  than  others  when 
the  mesh  is  coarse.  Examples  of  this  are  presented  later. 

One  catastrophe  that  can  occur  in  PIC  calculations  arises  when  the  cell 
next  to  a  left  reflective  boundary  has  an  empty  neighbor  on  its  right.  The 
energy  Equation  (8)  for  the  cell  becomes 


91  _  pu 
9t  ps 

and  the  acceleration  Equation  (7)  becomes 


9u  p 
9t  ~  ps 

If  u  and  I  (and  thus  p)  are  negative,  then  the  cell  is  unstable  since  both  9l/9t 
and  9u/8t  are  negative.  Starting  from  the  moment  of  zero  velocity,  the  time 
required  for  the  temperature  and  velocity  to  reach  — < «  is  T„  =  -f-  /  2s  \ 

2  V  (T_1)uo  j 

where  we  have  assumed  a  polytropic  gas  of  exponent  y  and  an  initial  accel¬ 
eration  in  the  cell  of  —  Uq  . 

Equation  (28)  shows  that  negative  temperatures  can  arise  if  the  cell 
size  is  large  enough.  Suppose  that  I  and  p  are  both  zero.  Then  9l/9t 
can  be  negative  if 
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91  s 
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9_ 

9x 


(pu  i) +  pu 

With  u  >  0  and  9l/9x  <  0,  this  condition  becomes 


<  0 


/  81  \ 

/ 9u\2  2 

91 

r  h)  "  - pu 

_  \  9x)  +  s 

8x 

9_ 

9x 

In  order  to  produce  negative  temperatures.  As  s  becomes  small  there  is 
reached  a  value,  Sq,  corresponding  to  a  particular  continuous  configuration, 
such  that  for  smaller  s,  dl/dt  Is  no  longer  negative.  Physically,  this  corre¬ 
sponds  to  passing  the  point  beyond  which  internal  energy  is  transported  into 
a  cell  faster  than  it  can  be  converted  to  kinetic  energy.  If  9l/9x  =  0,  then 
the  requirement  for  dl/dt  <  0  becomes,  instead, 


9  I  _  /  9u_ 

a  2  \9x 

9x 


V 


which  means  that  I  must  be  locally  a  maximum.  Since,  however,  1  =  0,  we 
see  that  this  situation  does  not  arise. 

A  further  condition  (more  stringent  than  necessary)  is  derived  by  setting 
9u/9x  =0.  In  that  case  we  see  that 


2p 

91 

9x 

At 

M  ) 

9x\ 

3 

9x  / 

which  means  that  to  avoid  negative  temperatures,  s  must  be  smaller  than  the 
distance  over  which  pl9l/9xl  changes  appreciably. 

These  ideas  indicate  some  of  the  causes  and  cures  for  negative  tem¬ 
peratures  and  the  inaccuracies  they  symbolize.  Peculiar  equations  of  state 
and  other  fictitious  situations  can  also  cause  negative  temperatures. 

The  time  interval,  fit,  between  cycles  has  already  been  restricted  by 
requiring  the  energy  discrepancy  of  Equation  (20)  to  be  sufficiently  small. 

It  is  also  restricted  by  the  condition  (cfit)/s  <  1  which  is  required  for 
stability  of  the  difference  equations.  Experience  shows  that  even  a  slight 
violation  of  this  condition  leads  to  catastrophe.  A  further  restriction  is  the 
requirement  that  <  1,  where  umax  is  the  maximum  fluid  speed  in 

the  system.  This  mechanical  restriction  arises  from  the  fact  that  no 
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allowance  is  made  for  a  particle  to  travel  more  than  a  cell  width  in  any  one 
cycle.  The  following  conditions  probably  make  the  last  restriction  even 
stricter:  The  success  of  PIC  depends,  clearly,  upon  there  being  a  statistical 
averaging  effect.  It  is  reasonable  to  suppose  that  this  will  be  enhanced  if, 
during  a  certain  elapsed  problem  time,  the  particles  experience  many  differ— 
r  ent  orientations  relative  to  the  mesh.  This  will  occur  if  (umax6t)/s  «  1  so 
;  that  many  steps  are  required  to  move  a  particle  across  each  cell. 

The  difficulties  concerning  statistical  fluctuations  and  their  proper 
averaging  in  PIC  raise,  perhaps,  the  most  difficult  question  to  discuss  in 
desired  generality.  There  are  situations  in  which  one  would  not  expect 
averaging  to  proceed  well.  An  example  is  the  case  of  a  plane  steady  shock 
! hitting  and  reflecting  from  a  rigid  wall.  Theoretically,  the  fluid  behind  the 
reflected  shock  comes  to  smooth  flat  profiles  of  temperature  and  density,  and 
the  material  velocity  is  zero.  Actually,  we  find  that  PIC  leaves  the  fluid  in 
a  perturbed  state  with  significant  fluctuations  about  the  constant  values.  These 
fluctuations  move,  roughly,  with  sound  speed.  The  particles  cross  boundaries 
only  infrequently  if  the  number  of  them  is  small,  so  that  the  damping  due  to 
repartitioning  is  slow.  Actually,  the  fluctuations  may  grow  owing  to  the  length 
of  persistence  of  cell-to-cell  particle  number  differences. 

It  is  found  that  the  most  satisfactory  results  are  obtained  from  PIC  for 
those  problems  which  have  a  general  mass  speed  that  is  cqmparable  to 
sound  speed  or  greater.  Systems  with  Targe  positive  accelerations  are  also 
well-treated. 

Difficulties  are  generally  encountered  in  problems  with  small  perturba¬ 
tions  in  a  fluid  which  is  otherwise  nearly  at  rest  and  in  problems  in  which 
the  fluid  is  supposed  to  decelerate  to  persistent  small  speeds.  Further,  in 
accelerating  fluids,  if  material  speed  remains  much  smaller  than  sound  speed 
for  a  time  long  enough  for  sound  to  travel  across  many  cells,  then  poor  re¬ 
sults  can  be  expected  in  the  form  of  large  fluctuations  about  the  proper 
values.  A  number  of  these  difficult  cases  are  illustrated  in  Chapter  m. 

The  origin  of  these  fluctuations  lies  in  the  instability  of  the  basic  ex¬ 
plicit  Eulerian  difference  equations.  As  the  fluid  speed  becomes  small  rela¬ 
tive  to  the  mesh,  Ad  of  Equation  (30)  decreases,  and  the  equations  approach 
the  unconditionally  instable  form.  Small  perturbations  tend  to  grow;  how¬ 
ever,  these  produce  local  velocities  which  in  turn  increase  the  effective 
viscosity,  preventing  unbounded  growth.  If  the  number  of  particles  per  cell 
is  increased,  the  perturbations  do  not  grow  so  large,  but  no  amount  of  in¬ 
crease  of  particle  number  can  cure  the  difficulties.  To  effect  a  complete 
cure,  we  need  a  stabilizing  mechanism  which  persists  to  zero  speed.  The 
requirement  is  for  a  dissipative  or  conductive  process  whose  cumulative 
effect  is  proportional  to  cell  size. 
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An  experiment  with  a  crude  approximation  to  heat  conduction  for  the 
most  pathological  cases  was  fairly  successful,  and  a  more  refined  technique 
of  using  this  or  similar  stabilizers  may  make  calculations  of  marginal  prob¬ 
lems  possible.  Such  artificial  procedures  are  not  required  for  the  large  class 
of  problems  to  which  the  method  is  suited.  None  of  these  smoothing  methods 
was  used  in  any  of  the  problems  reported  in  this  paper. 

3.  Several  Materials.  Another  problem  concerns  the  pressure  calcu¬ 
lations  when  two  or  more  materials  are  present  in  the  system.  Suppose  that 
a  cell  of  volume  r  (one-  or  more -dimensional)  has  particles  of  various  ma¬ 
terials,  a,  which  have  total  masses  Each  material  occupies  a  fraction, 

cra  (to  be  determined),  of  the  volume  of  the  cell. 

V  a  =  1  (32) 

Z_/  a 
a 


For  each  material  there  is  an  equation  of  state  of  the  form  p  =  f a(p,  I)  which 
can  be  written 


P 


a 


(33) 


Now,  in  reality,  across  any  material  discontinuity  the  pressure  is  continuous 
except,  possibly,  for  certain  discrete  moments  of  time.  This  leads  us  to 
postulate  that  the  pressures  of  all  the  materials  in  the  mixed  cell  are  always 
the  same.  Equating  these,  together  with  using  Equation  (32),  gives  just 
sufficient  information  to  find  the  unknown  quantities  cr^. 

As  an  example,  consider  all  the  materials  to  be  polytropic  gases  for 

which 


Pa  = 


(y  —  l)m  I 
a  a  a 


or  t 
a 


Equating  all  the  p^'s  to  p  gives 


cr  = 
a 


(y  -  l)m  I 
a  a  a 


P  t 


so  that,  from  Equation  (32), 
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(y  —  1)  m  I 
^  a _ a  a 

a  T 


In  this  case  the  total  pressure  is  just  equal  to  the  sum  of  the  partial  pres¬ 
sures  that  would  be  observed  if  the  gases  intermixed,  each  filling  the  whole 
cell.  This,  of  course,  would  be  true  for  any  system  of  materials  whose 
pressures  are  directly  proportional  to  their  densities. 

Returning  to  the  general  case,  we  now  specify  that  all  the  materials  in 
the  mixed  cell  have  the  same  velocity.  In  reality,  only  the  normal  compo¬ 
nent  of  the  velocity  is  continuous  across  a  material  discontinuity.  Forcing 
the  tangential  component  also  to  be  continuous  is  equivalent  to  introducing  a 
shear  viscosity,  which  will  be  discussed  in  more  detail  in  Chapter  n. 

In  Phase  I  the  rate  of  change,  9r/9t,  of  the  entire  internal  energy  of 
the  cell  is  calculated  using  an  equation  analogous  to  Equation  (8).  This  is  to 
be  distributed  among  the  materials  such  that 
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dr 

dt 


(34) 


(A  separate  temperature  is  kept  for  each  material  in  a  mixed  cell.) 

From  the  energy  conservation  equation  it  is  seen  that  if  pressures  and 
velocities  are  the  same  for  all  materials  in  a  close  vicinity,  then  so  is 
pdl/dt.  Thus,  with  A  independent  of  a, 
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=  A 


(35) 


Putting  this  into  Equation  (34)  we  obtain 

*  1  9T 

A  =  —  — 
t  9t 


or 
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(36) 
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Thus,  each  material  in  a  mixed  cell  gets  a  different  specific  Internal  energy, 
allowing  contact  temperature  discontinuities  to  arise. 

In  Phase  m  the  Incoming  particles  share  their  internal  energy  only 
with  others  of  the  same  kind.  In  the  velocity  adjustment,  each  material 
shares  in  the  internal  energy  increase  in  an  amount  proportional  to  its  final 
mass  in  the  cell. 
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Chapter  n 


EXTENSION  OF  METHOD  -  TWO  SPACE  DIMENSIONS 


Several  problems  arise  with  the  Introduction  of  a  new  degree  of  free¬ 
dom  for  the  fluid  motion.  Some  of  them  are  demonstrated  by  a  plane,  car¬ 
tesian,  two-dimensional  example.  Others  arise  only  with  the  introduction  of 
other  coordinate  systems. 


A.  The  Plane,  Cartesian,  Two-Dimensional  Box 


The  box  is  divided  into  a  system  of  square  cells  of  side  length  s,  in 
which  the  fluid  is  represented  by  particles.  We  assume  for  now  the  presence 
of  only  one  kind  of  fluid. 

Cell  number  ^  is  the  one  in  the 

lower  left  corner.  Cell  number  ^ 


has  center  coordinates 
0  +  1/2)  s' 

(i  +  l/2)s  *  By  this  convention, 
which  we  use  consistently,  the  cell 
comers  have  coordinates  which  are 
integer  multiples  of  s.  The  velocity 


in  cell  number  is 


All 


other  previous  notation  is  applicable 
here . 

The  density  in  cell  ( | )  is 


/ON 

VO/ 


(!) 
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(37) 


N jm 

2 

s 


so  that,  analogous  to  Equations  (7)  and  (8),  the  Phase  I  calculations  are 


Values  of 'uj,  'vj,  and  ij  are  obtained  as  in  Equation  (9). 

The  Phase  H  particle  motion  is  again  according  to  a  velocity  weighting 

procedure.  A  square  of  cell  size  is 
placed  around  each  particle,  and  the 
area  of  overlap  into  each  neighbor  cell 

_ determines  the  fraction  of  that  cell’s 

velocity  used  in  the  movement.  This 
is  done  for  each  component  of  velocity. 

I  I  Allowance  must  be  made  for  re- 

_ | . I - cording  the  passage  of  a  particle  into 

I  *  any  of  the  eight  neighbor  cells.  The 

1  I  repartitioning  of  Phase  ni  adjusts  the 

velocity  components  separately  in  order 

- - - that  no  change  in  vector  momentum 

occurs.  Thus,  for  example,  if  one 
particle  enters  cell  A  from  cell  B, 
then 
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The  associated  change  In  kinetic  energy  is 


<5KE  = 


N  A  +  1 

A 


rv  .  2  / — '  .  2 

<UA-V  +  VA  ~  VB 


so  that,  analogous  to  Equation  (13), 


NA  rA  +  ■b 
n“  +  1 

A 


1 6  KE| 
m(N^  +  1) 


Again,  if  particles  had  also  left  cell  A  in  this  cycle,  then  the  incoming 
particles  repartition  only  with  those  that  remain.  Clearly,  the  result  of 
repartitioning  for  several  entering  particles  is  independent  of  the  order  of 
treatment. 

The  one -dimensional  energy  considerations  of  Chapter  I  apply- to  this 
problem.  The  theoretical  discrepancy  is 


2  N3  m 
2  «  * 


Also,  boundary  conditions  are  determined  by  the  fact  that  each  border  cell 
has  an  energy  transfer  to  the  outside  of  the  system,  which  depends  upon  the 
pressure  and  normal  velocity  in  the  fictitious  outside  cell.  The  calculation 
is  analogous  to  that  leading  to  Equation  (21). 

The  treatment  of  empty  or  mixed  cells  in  this  example  follows  as  a 
direct  extension  of  the  procedures  of  Chapter  I. 
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With  regard  to  reflective  boundaries,  the  two-dimensional  problems 


introduce  the  following  new  difficulty. 


Suppose  that  a  rigid  obstacle  with  a 
comer  is  to  be  represented  by  a  re¬ 
flective  boundary  system  —  the  edge  of 
the  shaded  section.  With  fit  small 
enough  to  satisfy  the  accuracy  criteria, 
the  velocity  weighting  assures  that  no 
particle  will  ordinarily  enter  the  shaded 
section.  The  exception  arises  when  the 
particle  is  in  the  section  near  the  cor¬ 
ner  outlined  by  dots.  There  the  weight¬ 
ing  procedure  can  carry  the  particle 
into  the  shaded  area  no  matter  how 
small  fit  may  be.  To  overcome  this, 
a  process  may  be  used  whereby  the 

/  n+l\ 

y  i , 


new  coordinates  of  the  particle, 


/ 


(yn+1 

are  replaced  by  \  xn 

pending  upon  which  will  remove  the  particle  from  the  forbidden  zone. 

The  velocity  weighting  causes  an  effective  drag  area  to  be  formed  at 
slip  lines  between  fluids.  This  is  not  present  at  boundary  slip  lines  because 
of  reflectivity.  The  drag  is  not  dissipative  nor  does  it  propagate  into  the 
system  away  from  the  slip  line  as  long  as  the  slip  is  along  a  coordinate  line. 
For  diagonal  slips  it  is  not  possible  to  construct  a  straight  slip  line,  because 
at  least  some  characteristics  of  it  must  stair-step  along  coordinate  lines. 
What  is  really  calculated  in  this  case,  then,  is  approximation  to  the  flow  past 
a  rough  boundary.  This,  of  course,  will  send  a  signal  away  from  the  bound¬ 
ary  whether  or  not  viscosity  is  present. 

By  the  same  procedure  that  led  to  Equation  (27),  we  may  derive  the 
lowest  order  statistical  effect  of  PIC  in  this  plane,  2-D  system. 


P 


du 

dt 


=  —  Vp  + 


_9 

ax 


(46) 


where 


4 


The  form  of  this  result  shows  that  the  effect  of  repartitioning  is  not  directly 
analogous  to  either  shear  or  bulk  viscosity  in  the  usual  sense.  It  does, 
however,  have  some  of  the  properties  of  a  real  bulk  viscosity  but  none  of  a 
real  shear  viscosity.  There  is,  for  example,  dissipation  due  to  compression. 
Several  of  the  calculations  discussed  in  this  paper  demonstrate  the  effects 
of  repartitioning. 


B.  Problems  in  Other  Coordinate  Systems 


PIC  Is  adaptable  to  cylindrical,  spherical,  and  other  coordinate  systems, 
and  we  have  tested  a  variety  of  situations  in  which  the  configuration  is  inde¬ 
pendent  of  one  of  the  space  variables  (thus  2-D  problems).  In  these,  the 
particles  may  no  longer  be  mass  points,  but  more  complicated  geometrical 
figures,  each  having  a  different  fixed  mass.  A  number  of  difficulties  arise 
in  such  systems.  The  ways  in  which  one  can  write  space  differences,  de¬ 
termine  boundary  conditions,  effect  particle  movement  and  repartitioning,  etc . , 
have  been  tested  in  only  a  few  simple  cases.  The  manner  of  making  gener¬ 
alizations  is  an  extensive  subject  requiring  careful  treatment;  we  are  not 
presently  prepared  to  report  on  any  but  the  simplest  of  such  situations. 

The  use  of  cylindrical  coordinates  for  systems  with  azimuthal  symmetry 
*  will  illustrate  the  sort  of  difficulties  encountered.  Consider  in  particular  the 

energy  equation 


dE  I—.  .  -«►. 

p  dT“  -  v'(pu) 

(48) 

which,  in  cylindrical  coordinates,  may 
tially)  equivalent  forms 

be  written  in  either  of  the  (differed 

dE  1  3  (pur)  3(pv) 

^  dt  r  3r  3z 

(49) 

dE  _  pu  3(pu)  3(pv) 

^  dt  r  3r  3z 

(50) 

where  u  and  v  are  the  velocity  components  in  the  r  and  z  directions,  re¬ 
spectively.  When  these  are  written  in  difference  form,  they  are  no  longer 
equivalent.  A  criterion  for  choosing  between  them  is  based  on  the  considera¬ 
tion  of  a  spherical  system.  It  can  be  shown  that  one  form  of  difference 
equation  from  Equation  (49)  fails  to  preserve  sphericity,  the  discrepancy 
being  very  large  between  corresponding  cells  near  and  far  from  the  axis. 

The  same  form  from  Equation  (50),  however,  can  easily  be  shown  to  preserve 
sphericity.  (See  Chapter  IV,  Section  C.) 
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Boundary  conditions  can  be  derived  by  considering  the  energy  equation. 
With  cells  along  the  axis  labeled  j  and  those  In  the  radial  direction  labeled 
1,  the  change  in  energy  of  the  system  In  a  cycle  is,  using  difference  equations 
from  Equation  (50), 


6E  =  <5D  - 


2 

« 


tJ  St 


Hi  * 


2s 


(PU), 


(pv)J+1  -  (pv) 


j-1 

1 


(51) 


where  A  =  2irrj[  s2  (volume  of  cell).  The  axial  terms  cancel  In  pairs,  leaving 
boundary  conditions  at  the  ends  of  the  cylinder  which  are  analogous  to  those 
obtained  from  Equation  (21).  Employment  of  a  parts  summation  procedure 
for  the  radial  terms  shows  that  they,  too,  cancel  In  pairs;  and  in  terms  of 
a  fictitious  cell,  1  =  —1,  '*within»  the  axis,  conservation  of  energy  requires 
the  peculiar  condition 


(pu)3_l  =  <pu)j0 


(52) 


together  with  a  more  complicated  condition  at  the  outside  edge  of  the  cylinder. 

Arguments  involving  the  preservation  of  sphericity  show  that  the  effec¬ 
tive  velocity  for  moving  particles  must  be  derived  by  area-wise  weighting, 
as  in  the  plane  case,  with  no  weight  being  given  to  the  volume  of  rotation  of 
the  cell. 
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Chapter  m 


ONE-DIMENSIONAL  TESTS 


A  number  of  one -dimensional  tests  were  performed  to  reveal  the 
characteristics  and  limitations  of  PIC.  Most  of  the  difficulties  exposed  In 
these  tests  are  also  present  in  the  two-dimensional  applications.  The  most 
significant  results  are  summarized  In  this  chapter.  The  calculations  were 
performed  on  IBM  Electronic  Data  Processing  Machines,  types  701  and  704. 

A.  Simple  Steady-State  Shock 


PIC  was  coded  for  the  machine  calculation  of  a  shock  passing  through 
a  polytropic  gas  In  a  one -dimensional  box.  The  left  boundary  of  the  box 
allowed  for  the  Inflow  of  gas  in  the  shock  state.  The  further  progress  of  the 
shock  was  then  calculated  by  PIC.  Figures  1  and  2  show  the  appearance  of 
the  temperature  and  velocity  profiles  of  one  of  the  test  shocks  after  it  had 
traveled  a  number  of  cells.  Its  appearance  was  established  quickly,  and  then 
was  maintained  with  very  little  change.  There  were  four  particles  per  cell 
behind  the  shock  and  one  per  cell  ahead.  The  adiabatic  Index  was  y  =  5/3. 
The  curves  were  fitted  by  eye  to  pass  through  the  value  at  each  cell  center. 
The  shock  has  been  smeared  to  a  width  of  about  two  cells,  one  behind  and 
one  ahead.  (In  these  and  later  figures,  s  =  1.) 

As  a  more  stringent  test,  the  code  was  adapted  to  allow  calculation  for 
two  gases,  both  initially  cold  and  at  rest.  They  were  of  different  densities. 
The  same  shock  as  above  was  fed  in  from  the  left  and  allowed  to  develop 
its  steady-state  configuration  before  hitting  the  other  gas.  The  results  at 
late  times  are  shown  in  Figures  3  and  4  for  the  case  in  which  the  second 
gas  was  much  lighter  (rarefaction  reflected)  and  in  Figures  5  and  6  for  the 
case  in  which  the  second  gas  was  much  heavier  (shock  reflected). 

In  these  shock  problems,  the  smoothness  of  the  velocity  profiles  is  an 
indication  of  the  proper  entropy  production  across  the  shock.  The  smearing 
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of  the  shock  width  Is  no  greater  than  that  produced  by  other  finite  difference 
methods.  The  method  of  treating  contact  discontinuities  varied  in  slight  de¬ 
tail  from  that  which  is  described  in  Chapter  I. 

B.  The  Closed  One -Dimensional  Box 

A  machine  code  was  devised  which  would  calculate  the  configuration 
changes  of  one  or  two  fluids  in  a  closed  one-dimensional  box.  One  of  these 
fluids.  A,  was  a  polytropic  gas;  the  other,  B,  could  have  any  equation  of 
state  within  rather  wide  limitations. 

1.  The  Shock  Tube.  This  problem  constitutes  a  stringent  test  for  one¬ 
dimensional  calculation  procedures,  since  a  shock,  a  rarefaction,  and  a  con¬ 
tact  discontinuity  must  all  be  represented. 

Only  the  gas  was  put  in  the  tube;  initially  it  was  at  rest  and  uniform 
temperature.  A  diaphragm  separated  two  regions  which  were  at  different 
density  and  pressure  (we  take  the  high  density  to  be  on  the  left).  At  t  =  0, 
the  diaphragm  ceases  to  exist,  and  a  rarefaction  then  moves  leftward  while 
the  contact  discontinuity  and  a  shock  travel  to  the  right.  Beforejthe  dis¬ 
turbances  hit  the  walls  .of. jthe  box,  there  is  no  characteristic  length  in  the 
system,  so  thatlbe  effect  of  changing  cell  size  is  only  to  change  the  time 
scale.  The  only  artificial  parameter  that  can  be  varied  for  a  given  initial 
physical  situation  is  N,  the  average  number  of  particles  per  cell.  Several 
interesting  features  are  observed  from  the  resulting  velocity  profile: 

(a)  The  initial  perturbations  in  the  system,  resulting  from  the  way  In 
which  the  problem  was  started,  introduced  an  oscillatory  pattern  superimposed 
on  the  correct  velocity  profile.  This  means  insufficient  entropy  production 

in  the  early  stages.  The  fluctuations  travel  with  the  shock  front,  not  in¬ 
creasing  in  amplitude.  At  any  fixed  position,  they  die  out,  and  the  velocity 
comes  to  its  proper  value.  These  results  are  shown  in  Figure  7  for  the 
initial  configuration  of  a  5-to-l  density  ratio,  y  =  5/3. 

(b)  Increasing  the  number  of  particles  per  cell  everywhere  (maintaining 
the  same  initial  density  ratio)  reduces  the  amplitude  of  the  fluctuations.  In 
(a),  there  were  intially  five  particles  per  cell  on  the  left  and  one  per  cell  on 
the  right.  Multiplying  these  numbers  by  two  produces  roughly  one-half  the 
amplitude  of  oscillations  at  any  stage  of  the  calculation  but  otherwise  leaves 
the  rarefaction  appearance,  plateau  height,  and  shock  speed  the  same. 

(c)  Smaller  initial  density  ratios  across  the  diaphragm  produce  rela¬ 
tively  greater  oscillations.  A  problem  with  2-to-l  ratio  (two  particles  per 
cell  on  the  left,  one  per  cell  on  the  right)  produced  tremendous  oscillations, 
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200  to  300  per  cent  of  the  true  velocity.  Multiplying  the  number  of  particles 
by  six  reduced  the  fluctuations  to  an  amplitude  of  less  than  50  per  cent  of 
the  true  velocity.  In  both  cases,  however,  the  oscillations  were  centered 
about  the  true  plateau,  and  the  rarefaction  and  shock  speeds  were  well  pro¬ 
duced.  This  fact  is  significant  in  that  it  demonstrates  that  PIC  is  stable  to 
very  large  fluctuations  and  that  these  fluctuations,  even  when  extreme,  tend 
in  the  average  to  produce  the  gross  configuration  features. 


2.  Gravity  Problems.  Again,  only  the  gas  was  in  the  tube,  this  time 
at  constant  density  and  temperature  and  at  rest  everywhere.  The  box  was 
100  cells  long  and  closed  (reflective)  at  each  end.  An  acceleration  field  of 
magnitude  g  was  applied  to  the  right.  Until  an  element  of  material  feels  the 
signal  from  the  wall,  it  accelerates  uniformly,  and  the  velocity  profile  re¬ 
mains  perfectly  flat.  Theoretically,  a  shock  is  formed  at  a  distance  c2/2g 
back  of  the  wall  (c  is  the  material  sound  speed).  At  the  left  end,  a  rare¬ 
faction  is  produced.  As  discussed  in  Chapter  I,  the  material  which  comes  to 
rest  at  the  right  end  will  not  be  treated  well;  we  do,  indeed,  observe  resid¬ 
ual  fluctuations.  The  smoothness  of  the  rarefaction  zone,  however,  depends 
upon  the  magnitude  of  the  acceleration.  If  this  is  small,  the  slow -speed 
fluctuations  have  a  chance  to  become  appreciable;  if  it  is  large,  then  no 
fluctuations  appear.  Smoothness  seems  to  require  c2/sg  <  10.  This,  un¬ 
fortunately,  means  that  at  the  end  where  material  piles  up,  a  shock  theoret¬ 
ically  will  be  formed  at  a  distance  back  of  the  end  somewhat  less  than 
1/10  s.  These  results  put  severe  restrictions  on  doing  gravity  problems 
with  PIC.  No  cure  for  these  difficulties  is  effected  by  decreasing  6t,  the 
time  for  a  cycle. 

Two  typical  examples  are  shown  in  Figure  8.  Even  in  the  fluctuating 
example,  the  curve  is  tolerably  smooth  at  the  higher  velocities  (>half  of 
sound  speed),  showing  that  perturbations  to  the  flow  do  not  grow  in  rapidly 
moving  systems.  The  smoother  curve  was  allowed  to  develop  four  times 
longer  than  at  the  time  shown.  By  then,  a  large  cavity  had  formed  at  the 
left,  and  the  rarefaction  wave  was  still  smooth.  The  shocked  region  was  by 
then  quite  thick  and  poorly  represented.  Profiles  at  typical  later  times  are 
shown  in  Figure  9. 

Experiments  also  showed  that  a  small  perturbation  introduced  into  the 
middle  of  the  box  grows  if  c2/sg  is  much  larger  than  around  10,  and  other¬ 
wise  it  damps.  Increasing  the  number  of  particles  per  cell  does  not  cure 
the  situation,  either  here  or  in  the  rarefaction. 

The  case  of  infinite  acceleration  was  accomplished  by  running  the 
problem  of  a  gas  escaping  into  a  vacuum;  the  results  showed  good  agree¬ 
ment  with  the  theoretical  rarefaction  profile. 


3.  Coarse  Zoning  Test.  A  clear  demonstration  of  the  qualitative  dis¬ 
crepancies  that  can  occur  when  zoning  is  coarse  is  given  by  the  following 
example: 

The  box  is  in  three  sections;  from  left  to  right  are  fluid  A,  fluid  B, 
and  vacuum.  Initially,  the  two  fluids  have  the  same  pressure  and  both  are 
moving  to  the  right  with  the  same  speed.  One  expects  fluid  B  to  rarify  at 

first,  then  compress,  and  heat  as  it 
crowds  against  the  right  end.  A  shock 
A  I B  |  I  should  then  propagate  back  through 

~  |[  1  fluid  A.  These  qualitative  features  are 

observed  if  the  zoning  is  such  that 
fluid  B,  when  most  compressed  against 
the  end,  occupies  more  than  one  cell.  If,  however,  fluid  B  is  all  compressed 
into  one  cell,  then  a  vast  qualitative  difference  is  observed  from  that  ex¬ 
pected.  The  cell  of  fluid  A  next  to  the  hot,  crowded,  fluid-B  cell  drops  in 
temperature  when  it  should  rise;  its  temperature  becomes  negative.  No 
shock  pulse  is  sent  back  into  fluid  A  until  the  negative  temperature  is  suf¬ 
ficiently  raised  by  repartitioning  of  the  returning  fluid  B,  but  this  is  much 
too  late.  These  features  are  as  predicted  in  Chapter  I. 


C.  Expanding  Spheres 


Two  companion  spherically-symmetric  calculations  (one-dimensional  in 
spherical  coordinates)  were  run  under  the  following  code  names: 

OBOE,  using  the  spherical  form  of  PIC. 

ROBOE,  using  well  known  Lagrangian  hydrodynamic  methods. 

The  problems  started  with  the  following  configurations.  Three  concentric 
spheres  of  radii  5,  20,  and  30  marked  the  boundaries  separating,  respectively, 
a  central  piston  (not  zoned),  perfect  gas  A,  perfect  gas  B,  and  a  vacuum. 

The  gases  were  initially  cold  and  at  rest.  Gas  A  was  much  less  dense  than 
gas  B. 

The  situation  was  generated  by  allowing  the  piston  to  expand,  at  first 
accelerating,  then  with  nearly  constant  velocity,  then  decelerating. 

ROBOE  was  run  first;  some  of  the  observed  features  are  shown  in 
Figures  10  and  11.  In  order  that  the  OBOE  problem  have  the  same  input  at 
the  piston  the  following  procedure  was  adopted:  The  •'piston  cell"  was  de¬ 
fined  as  the  cell  with  the  piston,  except  that  if  that  cell  has  no  particles, 
then  the  piston  cell  included  also  the  next  cell  beyond.  The  input  data  in¬ 
cluded  both  the  piston  speed  and  the  total  energy  of  ROBOE  as  functions  of 
time.  The  particles  in  the  piston  cell  were  forced  to  have  the  piston  speed. 
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Knowing  the  kinetic  energy  input  that  this  produced  per  cycle,  we  subtracted 
this  from  the  total  energy  that  should  have  been  added  according  to  ROBOE 
and  Inserted  the  difference  as  Internal  energy  Increase  to  the  particles  In 
the  piston  cell.  This  was  equivalent  to  having  a  heat  conduction  at  the  OBOE 
piston  to  Insure  the  same  boundary  conditions  there  as  in  ROBOE.  This  did 
not  in  itself  force  agreement  between  OBOE  and  ROBOE  in  any  other  respects. 

The  energies  of  OBOE  are  shown  In  Figures  10  and  11,  together  with 
those  of  ROBOE.  Note  the  large  relative  fluctuation  of  internal  energy  in 
OBOE  at  early  times.  The  momentary  negative  internal  energy  arose  be¬ 
cause  the  true  thickness  of  material  which  should  have  the  piston  speed  was 
somewhat  less  than  the  width  of  the  piston  cell. 

There  is  some  difference  between  the  times  at  which  the  kinetic  energy 
in  gas  B  began  to  rise.  This  can  be  traced  directly  to  the  fact  that  the 
shock  front  was  steeper  in  OBOE  than  in  ROBOE.  The  Initial  rise  in  B 
kinetic  energy  in  ROBOE  was  due  to  the  precursor  ahead  of  the  shock.  This 
is  well  shown  by  Figure  12  where  the  velocity  profiles  for  the  two  problems 
are  shown  at  a  time  just  before  the  shock  hit  B,  and  by  Figure  13  which 
presents  the  motion  of  the  A-B  boundary  for  the  two  problems  as  a  function 
of  time.  It  is  not  known  which  calculation  more  closely  reproduced  the  true 
physical  picture. 

The  two  temperature  profiles  are  shown  in  Figure  14  for  a  time  which 
is  the  same  as  in  Figure  12.  Figures  15  and  16  show  velocity  and  tempera¬ 
ture  profiles,  respectively,  for  a  rather  late  time  in  the  problems. 
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Chapter  IV 


TWO-DIMENSIONAL  CALCULATIONS 


Various  two-dimensional  calculations  have  been  performed  with  PIC; 
several  of  those  which  demonstrate  the  characteristics  and  capabilities  of  the 
method  are  discussed  here.  An  unfortunate  feature  of  these  problems  is 
that  no  solutions  for  comparison  are  available  except  in  the  simplest  cases. 


A.  The  SUNBEAM  Code 

Two  different  polytropic  gases  are  arranged  in  arbitrary  configuration 
in  a  plane,  two-dimensional,  rectangular  box,  the  sides  of  which  are  rigid 
(reflective)  boundaries.  SUNBEAM  then  calculates  the  changes  of  configura¬ 
tion  in  the  manner  discussed  in  Chapter  n.  A  number  of  different  kinds  of 
calculations  can  be  done  with  such  a  code.  Two  examples  follow. 

1.  Interface  Motion  Studies.  Figure  17  shows  the  rectangular  box,  to- 
gether  with  the  interface  between  gases  at  various  times.  Initially,  the  gas 
at  the  right  was  of  low  density,  but  very  hot  so  as  to  have  a  high  pressure 
compared  with  the  fluid  to  the  left  which  was  of  greater  density  but  cold. 
There  were  four  particles  in  each  cell  at  the  start. 

A  number  of  variations  have  been  run.  No  apparent  difficulties  have 
arisen  in  the  calculations;  the  results  agree  with  the  qualitative  predictions 
that  one  can  make. 

2.  Cylindrical  Shock  Pulse.  Figures  18  to  21  show  a  sequence  of 
particle  configurations  resulting  from  the  following  initial  conditions:  The 
two  gases  (simulating  air  and  ground  except  that  the  density  ratio  was  1  to 
20)  were  both  at  the  same  pressure,  at  constant  densities,  and  cool  except 
that  within  a  circular  segment  of  five  cells  radius  in  the  upper  right  comer, 
the  temperature  was  100  times  greater  than  in  the  rest  of  the  air.  (Those 
cells  that  were  cut  by  the  circle  had  appropriate  fractions  of  the  large 
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central  temperature.)  Again,  there  were  no  computational  difficulties  in 
running  this  problem. 

B.  The  KAREN  Code 


A  single  polytropic  gas  is  confined  in  a  plane,  two-dimensional,  rec¬ 
tangular  box,  the  top  and  bottom  of  which  are  rigid  (reflective).  Horizontally, 
the  system  is  considered  to  be  one  period  of  an  infinite  periodic  system. 

This  is  accomplished  by  taking  for  the  conditions  in  the  fictitious  cells  just 
to  the  right  or  left  of  the  box,  the  conditions  just  on  the  Inside  of  the  corre¬ 
sponding  other  side.  Also,  particles  passing  out  across  one  boundary,  pass 
in  across  the  other.  Further,  a  fixed,  rigid  (reflective)  object  of  arbitrary 
shape  (as  long  as  its  boundaries  follow  cell  boundaries)  can  be  placed  any¬ 
where  in  the  box.  Initially,  the  gas  may  be  homogeneous  in  temperature, 
density,  and  velocity  (moving,  say,  to  the  right  at  a  specified  initial  Mach 
number),  or  the  gas  may  commence  motion  from  rest  under  a  horizontal 
acceleration  field.  A  large  number  of  problems  have  been  run  for  various 
^  initial  situations.  An  example  will  typify  the  kinds  of  results:  The  box  was 

square,  20  cells  on  a  side.  The  obstacle,  attached  to  the  bottom,  was  a 
rectangle  10  cells  high  by  6  cells  wide.  Initially,  the  gas  was  homogeneous, 

’  moving  to  the  right  at  near  Mach  2.  Figures  22  and  23  show,  respectively, 

the  horizontal  and  vertical  momenta  of  the  gas  as  functions  of  time. 

Problems  have  been  run  in  which  the  obstacle  was  placed  symmetrically 
in  the  center  of  the  box.  The  flow  remained  symmetric  until  a  slight  per¬ 
turbation  was  introduced,  whereafter  vertical  oscillations  in  momentum  grew 
to  much  greater  than  the  perturbation  intensity  and  showed  a  frequency  from 
which  a  very  reasonable  Strouhal  number  could  be  derived. 


C.  Cylindrical  Coordinates 


In  connection  with  the  cylindrical-coordinate  procedure  discussed  in 
Chapter  n,  we  mention  one  problem  that  was  designed  to  test  the  preserva¬ 
tion  of  sphericity.  This  was  the  calculation  of  an  expanding  sphere  of  hot 
gas  moving  under  its  internal  pressure  into  a  cold  gas  of  uniform  initial 
density.  After  the  central  gas  had  expanded  to  more  than  three  times  its 
original  radius,  the  deviations  from  sphericity  were  around  3  per  cent. 
Figure  24  shows  the  initial  boundary  between  the  gases,  together  with  the 
final  particle  configuration.  The  larger  dotted  line  Is  a  circular  segment; 
the  heavy  line  separates  particles  of  the  two  materials. 

With  this  coarseness  of  mesh  and  small  number  of  particles,  deviations 
are  expected  to  be  relatively  large.  Similar  calculations  with  much  finer 
zoning  and  many  more  particles  do  give  better  results. 
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Appendix  I 


OTHER  METHODS  FOR  SOLVING  2-D  PROBLEMS 


There  Is  a  large  difference  in  difficulty  between  1-D  and  2-D  hydro¬ 
dynamic  problems,  but  any  successful  procedure  for  solving  2-D  problems 
is  likely  to  give  relatively  little  difficulty  when  it  is  being  applied  to  the 
analogous  3-D  problems.  Here  we  discuss  briefly  some  other  methods  for 
solving  2-D  problems: 

1.  Moveable  Coordinate  (Lagrangian)  Methods.  The  essential  feature 
of  these  methods  is  that  the  coordinate  system  is  a  mesh  of  cells  imbedded 
In  the  fluid  and  moving  with  it.  A  typical  application  has  each  comer  of 
the  mesh  carrying  a  certain  fixed  mass  of  fluid  as  well  as  pressure,  position, 
and  velocity  which  vary  with  time.  The  appropriate  differential  equations  are 
written  in  a  finite  difference  approximation  form  and  advancements  of  the 
configuration  in  time  are  thereby  related  to  space  differences. 

The  Lagrangian  approach  has  proved  particularly  useful  for  treating 
systems  involving  several  fluids,  because  each  mesh  point  always  retains 
identity  with  its  initial  portion  of  the  fluid.  The  interfluid  boundaries  are 
always  clearly  delineated.  A  large  number  of  strikingly  successful  calcula¬ 
tions  have  been  performed  by  several  groups  of  workers. 

The  Lagrangian  methods  are  limited,  however,  to  use  with  systems  in 
which  no  large  distortions  of  the  fluid  occur.  In  a  finite-size  mesh,  various 
topological  catastrophes  can  happen  which  reduce  further  results  to  nonsense. 
Indeed,  rather  serious  doubt  is  cast  upon  the  accuracy  of  representing  the 
true  solution  when,  for  example,  a  system  whose  equations  were  based  on  an 
orthogonal  mesh  becomes  distorted  significantly  away  from  orthogonality. 

Problems  involving  oblique  collision  of  two  free  surfaces  are  difficult 
to  solve  by  Lagrangian  methods. 
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2.  Fixed  Coordinate  (Eulerian)  Methods,  The  essential  feature  of 
fixed  coordinate  methods  is  that  the  coordinate  system  is  not  tied  to  the 
fluid;  usually  it  is  stationary  in  the  laboratory  frame  of  reference.  Strict 
application  of  this  approach  allows  no  identification  of  the  fluid  elements. 

Each  cell  of  the  mesh  is  characterized  by  uniform  density,  pressure,  "color" 
(i.e.,  designation  of  material  kind),  and  velocity.  Thus,  finite  space  differ¬ 
encing  procedures  for  representing  the  differential  equations  can  retain  equal 
applicability  throughout  a  wide  range  of  fluid  distortions. 

The  principal  difficulty  with  the  strict  Eulerian  methods  is  that  they 
tend  to  introduce  false  diffusions,  especially  noticeable  with  material  boun¬ 
daries.  This  arises  from  the  fact  that  each  cell  is  forced  to  be  homogeneous; 
when  material  enters  a  cell,  its  characteristics  are  uniformly  mixed  with 
those  of  all  the  other  materials  in  the  cell. 

3.  Mixed  Euler-Lagrange  Systems.  Many  modifications  have  been  pro¬ 
posed— some  have  been  tried— which  take  advantage  of  the  better  features  of 
both  fixed  and  moveable  coordinates.  These  include:  (a)  Fixed  coordinates 

in  one  direction,  moveable  in  the  other.  Usually,  advantage  is  taken  of  the 
prior  knowledge  of  slip-line  positions,  (b)  Fixed  coordinates  but  with,  in 
addition,  a  special  Lagrangian  treatment  of  interfluid  boundaries,  (c)  Two 
superimposed  coordinate  systems— an  Eulerian  one  for  calculating  and  ad¬ 
vancing  field  functions  and  a  Lagrangian  one  for  keeping  track  of  the  positions 
of  material  points.  This  procedure  is  generally  somewhat  extravagant  in  the 
use  of  memory  space  in  the  computing  machine. 

4.  Repulsive -Particle  Methods.  The  fluid  is  represented  by  a  system 

of  mass  points  of  fixed  mass  whose  coordinates  vary  with  time.  The  pressure 
In  the  fluid  is  represented  by  prescribing  some  repulsive  interparticle  force 
which  is  a  function  of  separation. 

For  example,  the  pressure  in  a  polytropic  gas  in  n-dimensional  space 
can  be  represented  by  an  interparticle  repulsive  force 

_  —9  const. 

1)  =  8ru  [<yn<y-1) 

where  y  is  the  adiabatic  exponent  of  the  gas,  and  the  constant  varies  with 
entropy. 

Many  problems  arise  in  the  practical  application  of  these  and  any  other 
methods  for  solving  the  differential  equations  of  motion.  Some  of  the  diffi¬ 
culties  are  the  following. 
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1.  Stability  of  the  difference  equations.  Instabilities  may  arise  from 
several  sources  which  include:  (a)  Having  too  large  finite  time  interval  per 
cycle -the  classical  "Courant  instability. "  (b)  Failure  of  the  method  when 
an  attempt  is  made  to  treat  a  situation  of  a  type  which  the  method  is  not 
capable  of  handling,  (c)  Boundary  instabilities.  Free  surfaces  in  Lagrangian 
methods  may  be  troubled  by  instability.  The  method  presented  in  this  re¬ 
port  can  suffer  from  a  peculiar  boundary  instability  under  certain  conditions. 

2.  The  production  of  entropy.  In  order  to  allow  for  the  presence  of 
shocks  without  special  treatment,  a  dissipative  mechanism  must  be  present 
in  the  difference  equations.  This  is  usually  achieved  in  two  ways:  (a)  An 
artificial  viscosity  may  be  introduced  into  the  equations,  (b)  The  truncation 
errors  of  the  finite  difference  approximation  can,  under  certain  circumstances, 

automatically  provide  sufficient  dissipation. 

Some  success  has  been  achieved  by  actually  keeping  track  of  shock 
boundaries  separately  and  relating  field  functions  on  both  sides  by  the  proper 
conservation  conditions. 

3 #  Accuracy  of  representation  of  the  solution  of  the  differential  equations. 
There  is  no  guarantee  that  making  the  finite  difference  intervals  smaller  and 
smaller  will  produce  convergence  to  the  true  solution.  Many  authors  have 
devoted  a  large  amount  of  analysis  to  this  and  related  topics.  Every  slight 
variation  of  procedure  within  the  framework  of  any  one  method  must  be 
examined  for  its  contribution  to  inaccuracies. 
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Appendix  n 


AN  ALTERNATIVE  DERIVATION  OF  THE  PIC  DIFFERENCE  EQUATIONS 
by  Eleazer  Bromberg 


Introduction 

This  section  presents  a  derivation  of  PIC  which  emphasizes  the  role 
of  Lagrangian  and  Eulerian  concepts.  This  approach  may  help  to  develop 
variant  numerical  techniques  when  needed. 

We  shall  assume  for  simplicity  that  the  boundary  of  the  domain  D 
under  study  is  fixed  and  that  no  fluid  crosses  it;  the  resulting  equations  can 
easily  be  modified  to  take  account  of  different  boundary  conditions.  We  shall 
also  assume  that  we  are  dealing  with  a  single  material,  characterized  by  a 
single  equation  of  state  relating  the  values  of  density  p,  pressure  p,  and 
specific  internal  energy  I  of  the  material. 

As  in  the  customary  Lagrangian  formulation,  the  time  rate  of  change 
of  position  of  a  point  in  the  fluid  is  equal  to  the  velocity 

•  dx 

x  =  -=u(x,t)  (A-l) 


where  x  and  u  may  be  considered  as  scalars  in  one -dimensional  problems 
or  as  vectors  {xj}  and  {uj}  in  two  or  more  dimensions.  Integration  of 
Equation  (A-l)  yields 


x  =  x(x0,t)  (A-2) 


the  position  of  any  given  point  in  the  fluid  as  a  function  of  time  t  and  its 
initial  position  Xq. 
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A.  Mass  and  Density 


The  density  distribution  at  time  t,  p(x,t),  must  satisfy  the  mass- 
conservation  law  which  we  write  in  integral  form 

Jy  P(x,t)dV  *  X0PO(xO)dVO  <A-3) 

where  V  is  any  region  fixed  in  the  fluid,  V0  is  its  initial  volume,  p0  is  the 
initial  density  distribution,  and  dV  and  dVQ  are  corresponding  elements  of 
volume.  This  yields  the  alternate  form  of  the  mass -conservation  law 


P 


=  P 


% 
0  dV 


(A -4) 


where  dV0/dV  is  the  Jacobian  in  multidimensional  cases. 

We  now  introduce  the  assumption  that  p0  is  to  be  represented  by  a 
delta  function;  that  is,  by  a  function  which  differ®  from  zero  only  at  a 
finite  number  of  points  x^,  and  which  rises  to  infinity  at  each  of  those 
points  in  a  manner  such  that  the  volume  integral  of  the  density  in  any  suf¬ 
ficiently  small  neighborhood  enclosing  a  particular  point  is  a  constant 
characteristic  of  the  point  but  independent  of  the  size  or  shape  of  the  volume 
element.  This  constant  is  referred  to  as  the  "mass  of  the  particle"  as¬ 
sociated  with  the  point  in  question.  These  masses  are  chosen  so  as  to 
satisfy  the  conservation  law: 
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where  D  is  the  entire  domain,  which  is  considered  to  be  divided  up  into  n 
subdomains  such  that  each  subdomain  encloses  only  one  point  at  which 
pQ  is  other  than  zero.  For  any  arbitrary  volume  V,  we  obtain 


p(x,t)dV  = 


Vxo)dvo 


(A-6) 


where  the  sum  is  taken  over  those  particles  to  be  found  in  the  volume  V  at 
time  t.  This  volume  may  be  so  chosen  as  to  be  empty  of  particles  at  time 
t,  in  which  case  the  value  of  the  corresponding  integral  will  be  zero.  How¬ 
ever,  this  does  not  interfere  with  the  validity  of  the  conservation  law  in  the 
large . 


B.  Solving  the  Dynamical  Equations 

Having  established  the  characteristics  of  the  density  function,  we  can 
turn  to  the  momentum  and  energy  conservation  laws  and  develop  the  com¬ 
putational  process  based  on  the  above  considerations. 

We  assume  that  the  domain  D  is  covered  by  a  mesh  which  is  chosen 
so  as  to  lead  to  simple  expressions  for  the  difference  quotients  representing 
such  quantities  as  gradient,  divergence,  and  normal  derivative  associated 
with  any  cell  of  the  mesh.  The  mesh  distribution  does  not,  in  general,  have 
any  special  relation  to  the  subdomains  in  the  previous  section.  The  density 
associated  with  any  cell  c  is  obtained  from  Equation  (A-6)  by  approximating 
the  integral  on  the  left: 


pdV  =  pc 


V 


c 


from  which  it  follows  that 


2 


m 


(A-6*) 


(A-6**) 


In  the  process  of  integration,  we  treat  the  mesh  as  Lagrangian  at  the 
beginning  of  each  time  cycle,  allowing  each  cell  to  be  displaced,  but  then 
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we  immediately  rezone  so  as  to  return  to  the  original  spatial  mesh  distri¬ 
bution.  Thus,  each  time  cycle  of  integration  of  the  momentum  and  energy 
conservation  equations  is  divided  into  two  parts.  In  one,  we  deal  with  each 
cell  of  the  mesh  as  fixed  in  the  fluid  rather  than  in  space  to  determine  its 
change  of  momentum  and  energy;  this  displaces  its  boundaries.  In  the 
second  phase,  we  rezone  the  cells  (and  mass  points)  in  order  to  reintroduce 
the  original  mesh  distribution;  in  the  process,  we  repartition  momentum 
and  energy  so  as  to  obtain  new  spatial  distribution  functions  without  affecting 
the  conservation  of  these  quantities.  The  initial  mesh  structure  is,  there¬ 
fore,  kept  fixed  in  the  sense  that  it  reappears  at  the  end  of  each  integration 
cycle.  This  method  is,  consequently,  not  subject  to  the  difficulties  which 
normally  arise  in  the  Lagrangian  method  as  a  consequence  of  the  distortion 
of  the  Lagrangian  cells.  We  shall  henceforth  label  the  cells  considered  to 
be  fixed  in  the  fluid  Lj,  and  those  fixed  in  space  Ej. 

Phase  I 

The  momentum  and  energy  conservation  laws  can  be  written  as  follows: 


d_ 

dt 


d_ 

dt 


i 

X 


pu 


pe 


dV 


-i 


pn  dS 


dv  =  4  J  P<I  +  |u2)dV  =  -  J  (pu)n  dS 


(A -7m) 


(A-7e) 


where 


e  =  the  specific  total  energy 
I  =  the  specific  internal  energy 
dS  =  a  surface  element 
(pu)n  =  the  normal  component  of  pu 
n(x,t)  =  the  unit  outward  normal  to  the  surface 
SD  =  the  surface  bounding  volume  D. 

A  similar  set  of  equations  applies  to  each  cell  of  the  mesh: 


±1 


pu 


pn  dS 


(A-8m) 
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I  4/1  dV  +  I  •(,  7 p"2  dV  =  -  II,  (pU>n  dS  <A-8e> 

but  the  volumes  Li  and  their  surfaces  SLi  are  now  functions  of  time, 
though  they  coincide  with  the  fixed  mesh  at  time  t. 

Integrating  Equations  A-8  between  times  t  and  t  +  6t,  we  obtain 


f  pu  dV  =  f  pu  dV  -  r  f 

•'L^t+dt)  ^(t)  Jt  JS 


pn  dS  dt 


(A-9m) 


SLi(t) 


L 


L^t+fit) 


Idy  =  -- 


2  j*r 

pu  dV 


^^t+dt) 


-1 


2 

pu  dV 


Lt(t) 


J. 


+  I  pi  dV 

Lj(t) 


-CL 


SL^t) 


(pu)n  dS  dt  (A-9e) 


in  which  each  integrand  is  to  be  evaluated  over  the  indicated  domain  at  the 
indicated  time.  There  is  considerable  latitude  in  choosing  the  particular 
formula  for  evaluating  the  integrals  in  (A-9)  numerically.  In  the  present 
report,  we  use  the  following  prescriptions: 

(a)  The  integrals  over  time  are  set  equal  to  the  value  of  the  integrand 
at  the  time  corresponding  to  the  lower  limit  of  t,  multiplied  by  <5t;  that  is, 


pn  dS  dt 


(pu)  dS  dt 
n 


pn  dS 


(pu)  dS 
n 


(A-lOm) 


(A-lOe) 


(b)  The 


volume  integrals,  which  are  all  of  the  form 
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where  f  is  any  function  of  x  and  r  may  be  either  t  or  t  +  6t,  shall  be 
evaluated  as: 


pf  dV 


f  (r)  /  pdV 

L1  \(r) 


(A-ll) 


where  fj_,i(T)  is  the  value  of  f  associated  with  the  cell  Lj  at  time  r.  Since 
the  Lj  are  being  treated  as  Lagrangian  cells, 


the  summation  being  taken  over  the  mass  points  located  within  the  cell  L^. 

The  superscript  n  identifies  the  cycle  of  time  integration  (and  hence 
the  time  t).  We  use  the  following  nomenclature  for  f: 


(t  +  St)  =  f 
Li 


(A -12) 


The  numerical  integration  formulae  for  Equations  (A-8)  [or  (A-9)]  can, 
therefore,  be  written  as 


n  w>  n.  n  r 
M  u  =  M  u  —  St  /  pn  dS 

J  a 


(A-13m) 


n  ~  1  .  n 

M  I  =  --M 


~  2  n 
(u)  -  (u  ) 


+  MV-«  [ 

Jo. 


(pun).  dS 


(A-13e) 

The  subscript  Lj  is  xmderstood  to  be  associated  with  each  quantity  in 
(A -13) . 
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In  general,  the  values  of  the  Integrand  in  the  surface  integrals  must 
be  chosen  so  that  when  the  Equations  (A-13)  are  summed  over  all  the  cells 
Lj  of  the  domain  D,  the  sum  of  these  integrals  shall  yield  just  the  surface 
integral  over  the  bounding  surface  SD.  This  implies  that  the  contributions 
of  the  surface  integrals  shall  cancel  at  interfaces  of  adjacent  cells  (even 
if  one  or  more  of  these  cells  are  empty  of  mass  points).  In  the  present 
programs,  the  values  of  the  integrand  are  determined  from  the  relations 


[P] 


at  surface 


£  9p 
2  8n 


[pun] 

at  surface 


’  PL,  <uLl>n 


+  1  ±  (pu  ) 
2  8n  w  n' 


(A-13') 


Each  normal  derivative  is  represented  by  a  centered  difference  quotient, 
that  is,  by  the  difference  in  the  values  of  the  quantity  in  properly  chosen 
opposite  neighbors  divided  by  the  distance  between  the  centers  of  these 
cells.  Account  must  be  taken  of  the  fact  that  (uLi)n  ^ave  the  same 
absolute  value,  but  opposite  sign,  for  any  two  opposed  bounding  surfaces  of 
a  cell. 

Phase  n 

1.  Rezoning  of  Cells.  Having  determined  the  change  of  velocHy  and 
specific  internal  energy  associated  with  the  Lagrangian  cell  Lj  during  the 
interval  6t,  we  wish  to  rezone;  that  is,  to  determine  corresponding  quan¬ 
tities  associated  with  cells  (which  we  shall  call  Ej)  occupying  the  same 
spatial  positions  as  Lj(t)  [or  Lj(o),  for  that  matter]. 

For  this  purpose,  we  note  that  for  any  function  g(x,  t)  defined  over  a 
cell  Lp  at  any  time  t, 


d_ 

dt 


dS 


(A-14) 


where  E^  is  taken  to  be  a  fixed  volume,  equal  to  Lj(t),  and  SEj  is  the  sur 
face  of  Ej. 

Integrating  Equation  (A-14)  with  respect  to  time,  we  obtain 
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f  pg  dV  -  f  pg  dV  =  J  p(x,  t+fit)  g(x,  t+<5t)  dV 

J t  * t,  m  •'E. 


Li(t+6t) 


p(x,t)  g(x,  t)  dV 


t  •'SE. 


pg  dS  dt  (A-15) 


Since  Ej  =  L^t), 


pg  dV  =  /  p(x,t)  g(x,t)  dV 


(A -16) 


and  the  corresponding  terms  in  Equation  (A-15)  cancel. 

We  now  apply  Equation  (A-ll)  to  the  integral  over  Ej,  and  we  use  the 

labels 


gE  <t+st)  =  g^1 
i  1 


/  p(x,t+6t)  dV  =  2  m  =  M^+] 

E.  E  J  i 

l  i 


The  expression  (A-15)  can  be  written  as 


(A-17) 


n+1  n+1  H  rvj 

mEi  % 


pg  u  dS  dt 
n 


(A-18) 
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*• 


or,  since 


n+1  n+1  n 

M  S  —  M 
E.  SE.  E. 

11  i 


gT 


m  t+<5t  * 

J  t  *4e. 


pg  u  dS  dt 
n 


(A -19) 


The  surface  integral  corresponds  to  the  transport  of  g  from  one  cell 
to  the  other,  resulting  from  the  motion  of  the  fluid.  As  in  Phase  I,  it  must 
be  evaluated  in  a  manner  such  that  the  contributions  of  adjacent  cells  over 
their  common  interface  are  equal  and  opposite.  This  can  be  done  in  various 
ways;  in  the  current  computations,  we  replace  g  by  gL^  so  that  the  surface 
integral  becomes: 


(A -20) 


where  the  first  integral  on  the  right  hand  side  is  taken  over  those  parts  of 
the  surface  where  un  is  positive  (the  flow  is  out),  and  the  second  is  taken 
over  those  parts  where  the  flow  is  inward.  The  quantity  gjjjN  ^he 

neighboring  cell  from  which  the  flow  is  entering.  The  integrals  yield  the 
mass  of  material  entering  or  leaving  Ej  in  accordance  with  the  recipe  pre¬ 
scribed  in  Equation  (A-6).  This  can  be  seen  more  readily  by  replacing  ur  dt 
by  dxn;  the  latter  is  the  component  of  the  displacement  of  the  fluid  in  time 
dt  normal  to  the  surface  element  dS.  The  product  cb^  dS  is  the  volume  of 
fluid  passing  through  dS  in  time  idt.  A  typical  integral  can,  therefore,  be 
written  as 


ft 


(A-21) 
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where  6  V  is  used  to  indicate  the  special  character  of  this  element  of  volume. 
Transforming  in  accordance  with  Equation  (A-6),  we  get 


'>o6vo 


£ 

SE 


m. 


i,  fit 


(A-22) 


that  is,  the  sum  of  the  masses  of  the  particles  whose  location  has  crossed 
SEj  in  time  fit. 

Summarizing  Equations  (A-19),  (A-20),  and  (A -21),  we  obtain  the  typical 
rezoning  (or  repartitioning)  equation 


M. 


n+1  n+1 


n 


E. 

1 


g 


E. 

l 


=  M  g  -  g 
Ei  Li 


L. 

l 


£ 

out 


m. 


N 


g 


LjN 


£ 

in 


m 


j>  N 


(A -2  3) 


where  the  last  summation  is  over  those  neighboring  cells  from  which  particles 
have  entered  Ep  In  the  particular  cases  of  momentum  and  energy,  the 
corresponding  expressions  are 


M 


n+1 


n+1 

u 


M 


n+1 


n+1 

e 


(A-24m) 


(A-24e) 


The  appropriate  cell  subscripts  can  be  determined  from  Equation  (A-23). 
The  second  of  these  (A-24e)  is  generally  written  so  as  to  solve  for  In+1; 
namely, 


i 


l 
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2.  Relocation  of  Mass  Points.  In  order  to  determine  those  mass  points 
which  cross  any  cell  interface  in  time  fit,  it  is  necessary  to  keep  track  of 
the  positions  of  the  mass  points.  This  is  accomplished  by  integrating  Equa¬ 
tion  (A-l)  during  each  time  cycle: 


fiXk  =  \ 


St 


(A-25) 


or 


n+1 

*k 


=  Xk  + 


\ 


St 


(A-25’) 


where  the  subscript  k  is  used  to  identify  each  mass  point  and  as  its 
initial  position,  as  described  at  the  beginning  of  this  Appendix. 

It  is  important  to  note  that  the  Lagrangian  formulation  of  the  equations 
of  fluid  flow  (which  are  being  used  here)  has  as  one  of  its  basic  assumptions 
that  any  small  volume  element  fixed  in  the  fluid  remains  a  small  volume 
element  throughout  the  motion  and  always  includes  the  same  points  of  the 
fluid,  although  its  shape  and  size  may  be  altered  to  some  extent.  This  is 
equivalent  to  the  requirement  that  the  mass  points  do  not  pass  each  other  as 
originally  prescribed  in  the  description  of  PIC.  This  condition  is  satisfied 
if  fit  is  sufficiently  small  and  if  u(x,t)  is  represented  by  a  continuous  function 
when  used  to  determine  the  position  of  the  mass  points  as  in  Equation  (A-25) 
or  (A-25*).  The  value  of  u^  used  in  Equation  (A-251)  may  be  obtained  by 
some  form  of  interpolation  among  the  values  of  the  cell  velocity  for  groups 
of  neighboring  cells  using  either  uB,  ,  u  ,  u^+1,  or  some  combination  of 

±L»i  Ijj 

these  quantities.  An  effective  method  of  interpolation  is  presented  in  the 
body  of  this  report  (Chapter  n).  It  can  be  described  generally  by  the 
formula 
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(A-26) 


u 


“k  "L. 
o  x 


3 


3 

where  2  a  =  1,  the  cl  are  weight  functions,  and  L^j  refers  to  the  three 

*  «  k,  k, 

J=0  3  J 

neighboring  cells. 


C.  Summary 


The  computational  procedure  can  be  summarized  as  follows: 

We  start  with  a  suitably  chosen  mesh,  a  set  of  mass  points,  and  the 
values  of  the  mass,  velocity,  and  specific  internal  energy  associated  with 
each  cell.  (The  mass  associated  with  each  cell  must  be  equal  to  the  sum  of 
the  masses  of  the  mass  points  located  in  that  cell.)  Given  these  data  at 
the  end  of  the  n-th  time  cycle  (which  may  be  the  initial  time),  we  obtain  the 
corresponding  values  at  the  end  of  the  (n  +  l)st  time  cycle  by  solving  the 
following  systems  of  equations: 

1.  For  each  cell: 

P  =  P<I,p) 

the  equation  of  state,  where 

m  Mn 

o  =  V  J _ g- 

P  2->  y  y 

c  c  c 


(A -27) 


(A -6") 


where  V  is  the  volume  of  the  cell,  and 
c 


u 


n 

u 


6t 


n 


M 


pn  dS 


(A-13m) 


-  50  - 


(A-13e) 


in  which  the  surface  integrals  may  be  evaluated  by  any  one  of  several  methods 
described  elsewhere.  These  methods  must  always  be  consistent  with  over¬ 
all  conservation  laws  and  generally  also  have  to  provide  for  such  contin¬ 
gencies  as  "empty"  cells  and  special  types  of  boundaries  (such  as  the  axis 
of  symmetry  in  cylindrical  problems). 

2.  For  each  mass  point: 


n+1 

x 


n  ffx 

X  u  „  Ot 
eff 


(A -25") 


where  u  ..  has  the  form 
eff 


N 


"eff  *  “k  V  +  .2  “k.  UL  j 

o  i  3=1  J  i 


(A -28) 


N 

2  - 1 

J=0  j 


(A -2  9) 


and  u  is  a  velocity  associated  with  the  cell  in  which  the  point  is  located, 

LjO 

u  to  u.  _T  are  similar  velocities  associated  with  neighboring  cells,  and 


Lil  LiN 
the  a  are  weight  functions. 

K, 

J 


For  each  cell: 


Mn+1  *  Mn  -  2  mi  +  E  mi  N 
out  ^  in 


(A-30) 
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Velocity  profile  for  a  shock  having  hit  a  density  decrease.  Straight  lines  form  the  theoretical 
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Fig.  5  Temperature  profile  for  a  shock  having  hit  a  density  increase. 

Straight  lines  form  the  theoretical  profile.  The  two  curves  are  at 
different  times. 
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Fig.  7  Velocity  profiles  for  shock  tube  after  250  and  500  machine  cycles.  Position  and  magnitude  of 
peak  velocity  is  shown  for  everv  50  cycles.  Straight  lines  show  the  theoretical  profile. 


TIME 

Internal  energy  for  gases  A  and  B  for  the  problems  OBOE  and  ROBOE.  The  units  are 
arbitrary. 
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Fig.  12  Velocity  profile  in  OBOE  and  ROBOE  at  time  6-1/2. 
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Fig.  13  Position  of  the  A-B  boundary  as  a  function  of  time  for  the  problems  OBOE 
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Fig.  14  Temperature  profile  in  OBOE  and  ROBOE  at  time  6-1/2.  The  temperature  units  are 
arbitrary. 
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Fig.  15  Velocity  profile  in  OBOE  and  ROBOE  at  time  t  =  11-1/4.  B  and  B  are  the  A-B  boundary 
positions.  °  R 
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Temperature  profile  in  OBOE  and  ROBOE  at  time  11-1/4.  The  temperature  units  are  arbi 
trary.  Bn  and  BR  are  the  boundary  positions. 
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Fig.  17  Interface  positions  between  two  gases  at  various  times. 


Fig.  18  Particle  configuration  in  cylindrical  shock  problem  at  t  =  4.8. 
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Fig.  19  Particle  configuration  in  cylindrical  shock  problem  at  t  =  9.8. 
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Fig.  20  Particle  configuration  in  cylindrical  shock  problem  at  6  =  16.0 


Fig.  21  Particle  configuration  in  cylindrical  shock  problem  at  t  =  22.0. 
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Fig.  24  Particle  configuration  for  sphericity  test  problem  at  a  late  time. 


